// Filename: FFTDHTS.CPP

/************   Fast Fourier Transform   ********************\

1994 (c) CDI UI, The USSR Academy of Sciences, Samara
FFT function for double precision (symmetric).
    Uses huge array and table of exponents.

\************************************************************/

#include "fft.h"

int fftdhts( Comp16 huge *x, unsigned int Size, Comp16 huge *tab) {

  if (Size==1) return NoError;
  if (Size > FFTDHMAXSIZE || Size <= 0) return InvalidSize;

  int k = 1, Pow2 = 0;
  while ( k<Size ) { k *= 2; Pow2++; }
  if (k!=Size) return InvalidParams;

  for(int i=1; i<Size-1; i+=2) {x[i].re=-x[i].re; x[i].im=-x[i].im;}

  int j=0;
  for (i=0; i<Size-1; i++)   {
    Comp16 t;
    if (i<j) { t=x[j]; x[j]=x[i]; x[i]=t; }
    k=Size/2;
    while (k<=j) { j-=k; k/=2; }
    j+=k;
  }

  for (int l=0, ll=1; l<Pow2; l++, ll*=2) {
    Comp16 t,tv;
    for (j=0; j<ll; j++) {
      tv=tab[j<<(Pow2-1-l)];
      for (i=j; i<Size; i+=2*ll) {
	k=i+ll;
	t.re = x[k].re * tv.re - x[k].im * tv.im;
	t.im = x[k].re * tv.im + x[k].im * tv.re;
	x[k].re = x[i].re - t.re; x[k].im = x[i].im - t.im;
	x[i].re += t.re; x[i].im += t.im;
      }
    }
  }

  for(i=1; i<Size-1; i+=2) {x[i].re=-x[i].re; x[i].im=-x[i].im;}

  return NoError;
}
